11/15/2021

R packages

R-spatial




  • Open source! Free! (others include GRASS, QGIS)
  • Potentially more reproducibility relative to point-and-click workflows
  • Advantages to working within one program for your project from start to finish

Spatial data workflows

1. geospatial analysis and operations

Spatial data workflows

1. geospatial analysis and operations

2. visualizing spatial data

The {sf} package

  • integrates well with ‘tidy’ data workflows
  • vector-based spatial analysis made easy:

The {sf} package

  • integrates well with ‘tidy’ data workflows
  • vector-based spatial analysis made easy:
    • st_read() reads in spatial data (e.g., shapefiles)

The {sf} package

  • integrates well with ‘tidy’ data workflows
  • vector-based spatial analysis made easy:
    • st_read() reads in spatial data (e.g., shapefiles)
    • st_transform() transform or reproject spatial data

The {sf} package

  • integrates well with ‘tidy’ data workflows
  • vector-based spatial analysis made easy:
    • st_read() reads in spatial data (e.g., shapefiles)
    • st_transform() transform or reproject spatial data
    • st_buffer() buffer around data

The {sf} package

  • integrates well with ‘tidy’ data workflows
  • vector-based spatial analysis made easy:
    • st_read() reads in spatial data (e.g., shapefiles)
    • st_transform() transform or reproject spatial data
    • st_buffer() buffer around data
    • st_union() combine data into one layer

The {sf} package

  • integrates well with ‘tidy’ data workflows
  • vector-based spatial analysis made easy:
    • st_read() reads in spatial data (e.g., shapefiles)
    • st_transform() transform or reproject spatial data
    • st_buffer() buffer around data
    • st_union() combine data into one layer
    • st_intersection() crop one layer by another

Let’s take a closer look at some spatial data

Let’s take a closer look at some spatial data

Example field site: Swains Lake in Barrington, NH

Let’s take a closer look at some spatial data

Example field site: Swains Lake in Barrington, NH

Let’s take a closer look at some spatial data

Example field site: Swains Lake in Barrington, NH

Read in spatial data

Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):

#install.packages("sf")
library(sf)

local_lakes <- st_read("./data/NHDWaterbodyDataset_NHSeacoast.shp")

Read in spatial data

Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):

#install.packages("sf")
library(sf)

local_lakes <- st_read("./data/NHDWaterbodyDataset_NHSeacoast.shp")
class(local_lakes)

Read in spatial data

Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):

class(local_lakes)
## [1] "sf"         "data.frame"

Read in spatial data

Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):

head(local_lakes)

Read in spatial data

Load {sf} and read in a subset of the NH National Hydrography Dataset (available from www.granit.unh.edu):

head(local_lakes,3)
## Simple feature collection with 3 features and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -71.28925 ymin: 42.91414 xmax: -70.7096 ymax: 43.44024
## Geodetic CRS:  NAD83
##     Prmnn_I      FDate Resoltn GNIS_ID GNIS_Nm    AreSqKm Elevatn
## 1 141035025 2005-10-04       2    <NA>    <NA> 0.00100000      NA
## 2 141036768 2019-03-19       2    <NA>    <NA> 0.00028992      NA
## 3 141032014 2019-03-19       2    <NA>    <NA> 0.00038662      NA
##          ReachCd FType FCode VsbltyF      Shp_Lng      Shap_Ar
## 1 01060003011689   390 39004   24000 0.0010302955 7.489623e-08
## 2           <NA>   466 46600       0 0.0007273591 3.203682e-08
## 3 01060003008954   390 39004   50000 0.0007978068 4.298463e-08
##                         geometry
## 1 POLYGON ((-70.86617 42.9143...
## 2 POLYGON ((-71.28921 43.0608...
## 3 POLYGON ((-70.70971 43.4399...

Visually inspect spatial data

Interactive mapping tools are great for sharing and communicating spatial data, but also when you just want to look at the data.

One of these packages is {mapview}. Let’s take a look at our local_lakes dataset:

#install.packages("mapview")
library(mapview)

mapview(local_lakes)

Manipulate spatial data

{sf} integrates well with ‘tidy’ data workflows. What do I mean by that?

Manipulate spatial data

{sf} integrates well with ‘tidy’ data workflows. What do I mean by that?

Let’s subset our local_lakes dataset to only include my field site, Swains Lake:

library(dplyr)

swains_lake <- local_lakes %>%
  filter(GNIS_Nm == "Swains Lake")

Manipulate spatial data

{sf} integrates well with ‘tidy’ data workflows. What do I mean by that?

Let’s subset our local_lakes dataset to only include my field site, Swains Lake:

library(dplyr)

swains_lake <- local_lakes %>%
  filter(GNIS_Nm == "Swains Lake")

mapview(swains_lake)

Spatial Metadata

The st_read() function automatically stores information about our data, including:

  • the data format
  • the coordinate reference system (CRS)
  • the spatial extent of the data

Spatial Metadata

The st_read() function automatically stores information about our data, including:

  • the data format
  • the coordinate reference system (CRS)
  • the spatial extent of the data
st_geometry_type(swains_lake)
## [1] POLYGON
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

Spatial Metadata

The st_read() function automatically stores information about our data, including:

  • the data format
  • the coordinate reference system (CRS)
  • the spatial extent of the data
st_crs(swains_lake)
## Coordinate Reference System:
##   User input: NAD83 
##   wkt:
## GEOGCRS["NAD83",
##     DATUM["North American Datum 1983",
##         ELLIPSOID["GRS 1980",6378137,298.257222101,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4269]]

Spatial Metadata

Proper representation of coordinate reference systems (CRS) is very important for all spatial work!

CRS are often denoted using an “EPSG” code (see bottom of st_crs output). Common CRS include:

  • WGS84 (EPSG = 4326)
  • NAD83 (EPSG = 4269)

See other resources for more information, including https://datacarpentry.org/organization-geospatial/03-crs/

st_crs is your friend.

st_crs(swains_lake)

Spatial Metadata

The st_read() function automatically stores information about our data, including:

  • the data format
  • the coordinate reference system (CRS)
  • the spatial extent of the data
st_bbox(swains_lake)
##      xmin      ymin      xmax      ymax 
## -71.05927  43.18261 -71.02456  43.19972

Spatial Metadata

The st_read() function automatically stores information about our data, including:

swains_lake
## Simple feature collection with 1 feature and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -71.05927 ymin: 43.18261 xmax: -71.02456 ymax: 43.19972
## Geodetic CRS:  NAD83
##     Prmnn_I      FDate Resoltn  GNIS_ID     GNIS_Nm  AreSqKm Elevatn
## 1 141033653 2015-01-14       2 00870298 Swains Lake 1.379346      NA
##          ReachCd FType FCode VsbltyF   Shp_Lng      Shap_Ar
## 1 01060003002342   436 43600 1000000 0.1605085 0.0001527365
##                         geometry
## 1 POLYGON ((-71.05117 43.1887...

Spatial joins: find nearest weather station

Say we want to find the closest weather station to Swains Lake.

Spatial joins: find nearest weather station

Say we want to find the closest weather station to Swains Lake.

We can use the {rnoaa} package to identify NOAA climate stations, and then use spatial joins to select the nearest site and download some data.

library(rnoaa)

# Takes a minute to run:
stations <- rnoaa::ghcnd_stations() %>% 
  filter(state=="NH",element=="PRCP",last_year>2020) %>%
  st_as_sf(.,coords=c("longitude","latitude"),crs=4326) %>%
  arrange(id)

stations <- read_sf("./data/NH_GHCND_stations.shp") 

Spatial joins: find nearest weather station

Say we want to find the closest weather station to Swains Lake.

We can use the {rnoaa} package to identify NOAA climate stations, and then use spatial joins to select the nearest site and download some data.

head(stations,3)
## Simple feature collection with 3 features and 9 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -71.55 ymin: 43.4356 xmax: -71.3231 ymax: 43.5541
## Geodetic CRS:  WGS 84
## # A tibble: 3 × 10
##   id          elevation state name  gsn_flag wmo_id element first_year last_year
##   <chr>           <dbl> <chr> <chr> <chr>    <chr>  <chr>        <int>     <int>
## 1 US1NHBK0001      207. NH    TILT… <NA>     <NA>   PRCP          2009      2021
## 2 US1NHBK0002      214  NH    BELM… <NA>     <NA>   PRCP          2009      2021
## 3 US1NHBK0007      172. NH    LACO… <NA>     <NA>   PRCP          2009      2021
## # … with 1 more variable: geometry <POINT [°]>

Spatial joins: find nearest weather station

Calculate the distances between each weather station and Swains Lake.

st_distance(stations,swains_lake)

Spatial joins: find nearest weather station

Calculate the distances between each weather station and Swains Lake.

We need to transform the stations and swains_lake data so that they share the same CRS.

# Transform data to CONUS Albers Equal Area Conic projection:
swains_lake_proj <- st_transform(swains_lake,5070)
stations_proj <- st_transform(stations, 5070)

Spatial joins: find nearest weather station

Calculate the distances between each weather station and Swains Lake.

stations_proj2 <- stations_proj %>%
  mutate(dist = st_distance(.,swains_lake_proj,by_element = TRUE))

stations_proj2 %>% arrange(dist) %>% select(id,name,dist) %>% head(3)

Spatial joins: find nearest weather station

Spatial join: Find the closest weather station to Swains Lake.

closest_station <- st_join(swains_lake_proj,stations_proj,join=st_nearest_feature)

Spatial joins: find nearest weather station

Spatial join: Find the closest weather station to Swains Lake.

closest_station <- st_join(swains_lake_proj,stations_proj,join=st_nearest_feature)
head(closest_station)
## Simple feature collection with 1 feature and 22 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 1993730 ymin: 2506799 xmax: 1996588 ymax: 2508389
## Projected CRS: NAD83 / Conus Albers
##     Prmnn_I      FDate Resoltn  GNIS_ID     GNIS_Nm  AreSqKm Elevatn
## 1 141033653 2015-01-14       2 00870298 Swains Lake 1.379346      NA
##          ReachCd FType FCode VsbltyF   Shp_Lng      Shap_Ar          id
## 1 01060003002342   436 43600 1000000 0.1605085 0.0001527365 US1NHST0040
##   elevation state             name gsn_flag wmo_id element first_year last_year
## 1      80.8    NH BARRINGTON 3.2 E     <NA>   <NA>    PRCP       2016      2021
##                         geometry
## 1 POLYGON ((1994469 2507096, ...

Spatial joins: find nearest weather station

Now we can download and plot our precipitation data:

prcp_data <- rnoaa::ghcnd_search(closest_station$id, var = "PRCP")

# Plot data:
prcp_data$prcp %>% filter(date>"2021-08-01" & date<"2021-10-01") %>% 
  mutate(precip_mm = prcp/10) %>%
  ggplot() + geom_line(aes(x=date,y=precip_mm)) +
  labs(x = "Date (2021)", y = "Precipitation (mm)") + theme_classic()

Mapping vector and raster data

Useful packages for mapping in R:

  • {ggplot2} along with geom_sf()
  • {tmap} for customized, publication-quality maps
  • {mapview} and {leaflet} for interactive maps!
  • {rasterVis} for raster data

Mapping vector and raster data

We can still use {ggplot2} for plotting spatial data.

ggplot() + 
  geom_sf(data=swains_lake) + 
  coord_sf()

Mapping vector and raster data

We can still use {ggplot2} for plotting spatial data.

A lot of the formatting should already be familiar if you’re a {ggplot2} user:

ggplot() + 
  geom_sf(data=swains_lake,fill="cyan3",color="cyan4") + 
  coord_sf(ylim = c(43.18, 43.2)) + 
  ggtitle("Swains Lake, NH") + 
  theme_classic() 

Mapping vector and raster data

Let’s load one more spatial dataset to include in our map - a shapefile of the watershed boundary, that is, the area of land that drains into Swains Lake.

swains_watershed <- read_sf("./data/Swains_watershed_boundary.shp")
st_crs(swains_watershed)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

Mapping vector and raster data

Let’s load one more spatial dataset to include in our map - a shapefile of the watershed boundary, that is, the area of land that drains into Swains Lake.

ggplot() + 
  geom_sf(data=swains_lake,fill="cyan3",color="cyan4") + 
  geom_sf(data=swains_watershed, fill=NA, color="gray20", lwd=0.6) +
  coord_sf() + 
  ggtitle("Swains Lake, NH") + 
  theme_classic() 

Mapping vector and raster data

Add scales, orientation arrows, etc. using {ggspatial}.

You can probably think of other ways you want to edit your map for a report or publication, but we’re definitely getting closer!

library(ggspatial)

swains_map <- ggplot() + 
  annotation_map_tile(zoom = 15) +
  geom_sf(data=swains_lake,fill="cyan3",color="cyan4") + 
  geom_sf(data=swains_watershed, fill=NA, color="gray20", lwd=0.6) +
  ggtitle("Swains Lake, NH") + 
  theme_classic() +
  labs(x="Longitude",y="Latitude") + 
  coord_sf() +
  annotation_north_arrow(location="bl", 
                         height = unit(.75, "cm"),
                         width = unit(.5, "cm")) +
  annotation_scale(location="bl", 
                   pad_x = unit(1, "cm"), 
                   text_cex = 0.5, 
                   height = unit(0.2, "cm")) 

Mapping vector and raster data

Add scales, orientation arrows, etc. using {ggspatial}.

You can probably think of other ways you want to edit your map for a report or publication, but we’re definitely getting closer!